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The incremental unknowns-a multilevel scheme 
for the simulation of turbulent channel flows 

By M. Chen 1 , H. Choi 2 , T. Dubois 3 , J. Shen 1 AND R. Temam 4 


In numerical simulation of complex flows, it is important to identify different 
length scales of the flow and treat them differently. In this report, we introduce a 
new multilevel scheme for simulating turbulent channel flows. Two different versions 
of the scheme, namely the spectral and finite difference versions, are presented. The 
spectral version of the scheme is based on a spectral- Galerkin formulation which 
provides a natural decomposition of the flow into small and large wavelength parts, 
and which leads to linear systems that can be solved with quasi-optimal computa- 
tional complexity. In the finite difference version, the “Incremental Unknown” (IU) 
is used to separate the length scales. Preliminary numerical results indicate that the 
scheme is well suited for turbulence computations and provides results which are 
comparable to that by Direct Numerical Simulation (DNS) but with significantly 
less CPU time. 


1. Motivation 

The numerical simulation of turbulent flows is an extremely challenging task for 
both the numerical analysts and computational fluid dynamicists. The computing 
power required to resolve the enormous number of degrees of freedom and their 
nonlinear interactions involved in a turbulent flow is often near or beyond reach 
of the current computer capacity so that conventional numerical schemes are often 
impractical for turbulence simulations. 

The aim of this paper is to introduce a new multilevel scheme which is based on 
a differentiated treatment for small and large wavelength parts. It is well known 
in turbulence theory that the large number of small wavelengths only carry a small 
part of the total kinetic energy of the flow, however, the effect of their nonlinear 
interactions with large wavelengths over a long term integration can not be ne- 
glected and must be adequately resolved. Nevertheless, the small wavelength part, 
especially their nonlinear interactions, do not need to be represented in the same 
accuracy as the large wavelength part. Our multilevel scheme is specially designed 
such that it would produce results comparable to that by DNS but at significantly 
less cost so that one can simulate more complicated flows with limited capability of 
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the computer. The method can be applied to a class of dissipative equations and 
can be combined with a large number of existing numerical methods. 

The method starts with separating the length scales of the solution u as 

u - f -f g + r 

where / is the large length scale, g is the intermediate length scale, and r is the 
small scale. Then the different scales of the solution are treated differently, which 
could involve (a) neglecting some higher-order terms involving the small scales, (b) 
updating small scales with lager time interval. The effect of these further approxi- 
mations would, if done correctly, reduce the CPU time for each time step, improve 
the stability (the CFL condition will be only related to the large wavelengths), and 
allow larger time steps. 

There are two ways to look at this method. One is that we neglect some effect 
of the small scale terms. Another way is we think that large scale approximation 
is not enough, so we take into account the effect of small scale terms in an efficient 
way instead of simply adding more mesh points. 

This method has been applied to the simulation of 2D and 3D forced homogeneous 
turbulence (see Dubois, Jauberteau k Temam 1995a, 1995b, 1996 and the references 
therein). In the 3D case, it has been shown that the main statistical properties of 
homogeneous turbulence is well predicted with multilevel schemes. Indeed, while a 
saving in CPU time of 50-75% versus a classical Galerkin method is obtained, the 
energy and enstrophy spectra as well as the high-order moments of the velocity and 
its derivatives are accurately computed. The comparison of these results has been 
done with the results of direct simulations. 

In the case of homogeneous turbulence, when Fourier expansion of the velocity 
is used, the separation of the flow into large and small scales is trivial. However, 
this is not obvious for the channel flow problem because of the no-slip boundary 
conditions at the walls. In particular, the popular spectral-tau (Gottlieb k Orszag 
197/) method is not suitable for this purpose. We shall use the spectral-Galerkin 
method developed by Shen (1994, 1995) for the non-homogeneous direction. This 
spectral-Galerkin formulation not only provides a natural decomposition of the flow 
into small and large wavelength parts, but also leads to linear systems that can be 
solved with quasi-optimal computational complexity. 

In the finite difference case, we will use the IU’s developed by Chen k Temam 
(1991). The IU method has been used for steady equations, and the result is similar 
to preconditioning the associated matrix. The scheme was shown theoretically 
convergent and has an improved efficiency (Chen k Temam 1993). Here for the 
first time, the IU method is applied to unsteady problems. 

This report is an interim report: more detailed results using the new scheme for 
the turbulent channel flows will be reported later. 
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2* Incremental unknowns in the spectral case 

2.1 Formulation of the equations 
We consider the Navier- Stokes equations 

^ - uAu + (u ■ V)u 4- -VP = 0, (2.1) 

at p 

div u — 0, (2.2) 

in a channel ft = (0,L X ) x (-1,1) x (0 ,L Z ) with the boundary conditions: u = 
(u,u,u>) is periodic in x and z, and no slip on the walls. For this channel flow, we 
assume that the pressure P takes the form P = P + A> x, where P is periodic in 
directions x and z and Kp is a given constant. 

Following Kim, Moin &: Moser (1987), we set 


A — (ia * V)u — (A x ,Ay, A z ), 
du dw 

f = lh + Ik' 

du dw 
dz dx ’ 


d 

/ d\ x 

d\y\ 
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V dy 
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dz \ 

K dz 
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f dA z 

dAA 

= -{u 

•V) 9 

dv 

dv 

dw 

-( 

K dz 

dx J ‘ 

+ dy 

+ dx 

dy 


dv du 
Yz v 


(2.3) 


then, (2.1)-(2.2) are equivalent to the following equations (cf. Kim et ai 1987): 


— At? — uA 2 v — /i v (u,tt), 
dt 


dg 

dt 


- vAg = h g (u,u ), 


r dv 

/+ % =a 


(2.4) 


From the boundary conditions of u and the continuity equation (2.2), we deduce 
boundary conditions for v and g: 


^ il j 2 j / ) 0, 

g(x,±l,z,t) = 0. 

We emphasize that h„(-, •) and h g (-,-) are indeed bilinear forms since they are 
derived from the original bilinear form by linear differential operations. 
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Writing the Fourier expansion in directions x and 2 for u 


u(x,t)= Y, * = (*„*,), 

kez* 

where it j. = u’fc), and similarly for /, g and h v , h g , we derive from (2.4) 

that 

d,.; d 2 


*<* - 


dy 


j /2 


’fc + r/ (a 4 - 2A 2 |^ + j£f) 


dOj, 

5 fe(± i) =-^(± i )=°> 


(2.5) 


and 


d h 

dt 


d 2 

+ ^ k2 ~^h = h g, *(«’«)’ 


( 2 . 6 ) 


fffe(±l) = 0, 

where A 2 = (f^) A- 2 + (f^) A 2 . 

From the equations relating the velocity components u and w to / and g in (2.3), 
we derive 

27 r ^ 27T A ? 

ik x j~u k +ik z — W k = /jfc, 

2 ; 2 ; for all (*„*,) *(0,0). (2.7) 

***I7“ fc ” ikt j~™k - h’ 


For (fc x ,fc*) * (0,0), the relations (2.7) can be used to determine t/j,(y,J) and 


u^(y,t) in terms of f k {y,t) and g k (y,t). Hence, to complete the system, we still 
need additional relations for u 0 (y, t) and w 0 (y, t). To this end, we integrate the first 
and last components of the Navier- Stokes equations with respect to x and 2 to 
obtain 


duo d 2 iio 

~dt ~ V Ik? 



du 

v(x) — (x)dz + Kp = 0, 
dy 


du> 0 

~dT 


3 2 u> 0 

l/ w 


+ 



/ v dw . , _ 

v(x)—(x)dz = 0. 
dy 


( 2 . 8 ) 


The time discretization of (2.5), (2.6), and (2.8) is achieved by using a semi-implicit 
scheme with the second-order Crank-Nicolson for the linear terms and a third order 
explicit Runge-Kutta scheme for the nonlinear terms. Hence, we only have to solve 
a sequence of one-dimensional second-order equations for g k {y,t) and fourth-order 
equations for v k (y,t). 

Kim, Moin k Moser (1987) applied a Chebyshev-tau approximation to the 
y— direction. Since the direct application of tau method to fourth-order equations 
is unstable (Gottlieb k Orszag, 1977), they proposed a time splitting scheme which 
consists of solving several successive second-order problems to enforce the boundary 
conditions on v by using a technique similar to the influence matrix method. 

Based on a sequence of recent work by Shen (1994, 1995, 1996), we present below a 
spectral-Galerkin scheme for these second-order and fourth-order equations. Using 
this method, the system (2.5)-(2.6) can be directly solved. 
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2.2 A spectral- Galerkin approximation of the Kim- Moin- Moser formulation 
A Fourier- Galerkin approximation in the x and z directions is first applied to the 
problems (2.5) and (2.6), i.e. we look for 

t i N (x,t)= ^2 Uf l (y,t)e^ kx ^ x+kz ' I ^‘\ (2.9) 

(where N = (N x , N z ) and S/v = {k € Z 2 /(k x ,k 2 ) G [1 — x [1 — ^>S ^"]}) 

as a solution of the system of 


dt 




q2 / ^4 \ 

- w )ils + "V + w) 6k = 


*t(±il = ^ L (± 1 ) = o, 


( 2 . 10 ) 


and ^ 2 

w + ' ,{kt ~h )ik "' h ’- k{u, ‘' UN) ' (2.ii) 

$fc(±l) = °i 

for all k G S;v* 

We now describe Galerkin approximations of (2.10) and (2.11) in the y— direction. 
Let us denote 

• Pm : the space of polynomials of degree less than or equal to M, 

• Vm = span{v?(j/) G Pm ■ ^(±1) = 0}, 

• W M = span{<^(y) G Pm ■ vK±l) = 0, §*(±1) = 0}. 

Let pj(y) be either the Legendre or Chebyshev polynomial of degree j, then 


Vm = span{<£o, <t>\, • • • , <f>M- 2 } 


with <f>j{y) = Pj(y ) - p J+ 2 {y)- Moreover, following Shen (1996), we can determine 
( a j ibj) suc b that 

V'i(y) = pj(y ) + ajpj+2(y) + f>jPj+i(y) 

satisfies the boundary conditions t/’ j ( ± 1 ) = ^-(±1) = 0, i.e. V’j G W\j. Therefore 

W M = span{V>o,V’i,---,0itf-4}- 


The spectral- Galerkin scheme in the y— direction for (2.10) and (2.11) is to find 
vn,m(xJ) such that M (y,t) G Wm, and un,m( x ^) (similarly for w and g) such 
that Ufc M (y,t) G Vm , for all k G S/v, such that 




( [k 2 


d 2 

dy 2 




+ V 



d 2 (P 

dy 2 5t/ 4 




v J i {UN,M, u N,M 



U) 


( 2 . 12 ) 


for all j = 0, . . . , M — 4, 
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and 


(9 / ^2 \ 

^ + " ^( fc2 - Qyi'fik'ti) =(hgM UN - ^ 213) 

for all j = 0, . . . , M - 2, 


where ( 7 ?, t/Ow — /_!j ^(y ) 0(y) w dy with u>(y) = 1 in the Legendre case and uj(y) = 
(1 — y 2 )“? in the Chebyshev case. 

It is easy to see that in (2.13) the mass matrix M with entries m ; / = ( 0 /,^)^ 
is a sparse symmetric matrix with three nonzero diagonals, and that the stiffness 

o2 

matrix S with entries s 3 i = {-§^z is diagonal in the Legendre case, and is a 
special upper triangular matrix in the Chebyshev case such that the linear system 
(a Ad + S)x = b associated with (2.13) can be solved in 0(M ) operations (Shen 
1995). Similarly, the linear systems in (2.12) can be solved in O(M) operations, 
see Shen (1994, 1995). We emphasize that the above spectral- Galerkin scheme is 
superior, in both efficiency and accuracy, to the tau-method used in Kim, Moin & 
Moser (1987), and is, in particular, suitable for multilevel decomposition. 

The Legendre- Galerkin method has been implemented and tested. In this code, 
the pseudo-spectral computation of the nonlinear terms is done at the Chebyshev- 
Gauss-Lobatto points in the normal direction (see Shen 1996). A 128 x 129 x 128 
simulation at the Reynolds number of 180 has been conducted. The statistics have 
been compared to the one presented by Kim, Moin &; Moser (1987). 

2.8 A multilevel spectral- Galerkin scheme 

We now describe a multilevel scheme for the time integration of (2.12) and (2.13). 
For the sake of simplicity, we will only present a scheme based on a first-order semi- 
implicit scheme for the time discretization. However, one can easily generalize it to 
higher-order semi-implicit scheme. 

The basic idea of the multilevel scheme is to decompose the solution into several 
length scales and treat them differently in order to improve the efficiency and sta- 
bility of the classical Galerkin approximation. The special basis functions 
provide a natural decomposition of small and large wavelengths for this purpose. 
Furthermore, the small and large wavelengths are quasi- orthogonal in the following 
sense: 


{<Pl,<t>j) u, = o, for 


( d 2 4>i 
[ dy 2 ’ 
,d 2 h 
1 dy* ’ 


= 0, 

dj )w = 0, 


2 , 

for / 7 ^ j (Legendre case), 

for / < j or / + j odd ( Chebyshev case 


), 


(2.14) 
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H’lA’j L = 0, for j ± /, / ± 2, / ± 4, 

Fp-xU, 

(’dyT '^ “> = °’ for / f j,l ± j ± 2 , 

’ib I 

( yj-, 0>)u; = o, for l p j (Legendre case), 


(2.15) 


i 

( — = 0, for / < j or / + j odd (Chebyshev case). 

oy 

Given two appropriate cut-off numbers M p , M q such that 0 < M p < M q < A/, 
we may decompose u ^ jW (y,0 G Vm as follows 

M —2 

« fc,M (W'*) = 53 

i=o 

Mp-2 M,-2 M —2 (2.16) 

= “Jb.j^i(y) + 53 Ukj<f>j(y) + u k j( f) ] (y) 

J—O j= Mp - 1 j=M q — 1 

= Pu(y,0 + g«(y,0 + r u (y,0> 

and similarly for u>j, M (y,<) and then for g ^ M (y,t), for all fc G Sjv- Note that for 
the sake of simplicity, the dependence of p u , q u , and r u in k is omitted. We may 
also decompose M (y,f) € Wm as 

M- 4 

Sfc.Mfo’O = 13 '’k.j'J’Av) 

J=0 

M p —4 M q - 4 M-4 (2.17) 

= 53 + 53 5 4j^(w)+ 53 {, k,^j(y) 

j—Q j — Mp — 3 3 

= Pt-(y, 0 + 9i>(y, <) + r„(y, t). 

We finally obtain the following decomposition for Ufa M : 

Uk,M=P + < l + r ' 

where p = ( Pu,Pv,Pw ) and similarly for q and r. The decomposition (2.16) on u ^ M 
and Wfa M induces a decomposition of g ^ M into 

&,*(»•*) = Pff +99 + r 9- 

Then, thanks to (2.15) (resp. (2.14)), we can approximate the system (2.12) (resp. 
(2.13)) in Wm p (resp. Vm,) as follows 

| (V - + - ((*• “ 2t ’|> + |r) *•*). 

= (^,fc(p + q + r,P + q + r),^>) w , 

for j = 0, . . . , M p - 4, 
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Pg + ^' 2 (Pf,^i)u, - v 

= (^g,jk(P + 9 + » , -P + 9 + »’),^)^, 

for j = 0, . . . , M p - 2, 

and in W\f (resp. Vm ) as follows 


( U ' 2 “ ^ ){Pv + «•)’*>) + " ((** " 2fc2 |^ 


+ 5y 4 ) ^ Pu + 


= (*v,fc(p + 9 + »%P + 9 + 

for j = 0, . . . , M q - 4, 

^(p 9 + 9».*i) w + *>* (p» +g»,^>)« - K^t(p 9 + </ 9 ),<Aj)u, 


(2.20) 


= (j»,Jfe(P + 9 + r,p + q + r),<f>j) u , ^ 2-21 ^ 

for j = 0, . . . , M q — 2. 

Note that in (2.18)-(2.19) and (2.20)-(2.21) linear interaction terms coming from 
(Pg i )w (resp. (q v ,<f>j) u>) and similarly for r g (resp. r v ) are neglected. Until 
numerical tests are performed, it is not clear whether or not these terms have to be 
neglected. However, for the sake of simplicity we do not take them into account in 
the large or intermediate scale equations. 

By projecting (2.12) (resp. (2.13)) onto the space W M \W M% (resp. V M \V Mq ) we 
obtain the small scale equation 


dt\ {k ~d?^0 u +v 
= (Ufe(P + 9’P + 9)’V>>) I 




■ I (<* 2 - - v ( <l ' 4 - 2t ' 2 


( 2 . 22 ) 


for j = M q - 3, . . . , M - 4, 


d d 2 r 

fr( r 9’<l>ih + vk 2 {r g ,4>j) ' U - ^i-QyT^A 


(2.23) 


= ( A , t *(P + 9.P + 9).^)w- ft(q g ,<h)u,-i'k 2 (q g ,<t>j)„, 

for j = M q - 1, . . . , M - 2. 

We note that in (2.22)-(2.23) the nonlinear interaction between the small wavelength 
part t and the larger wavelength parts (p + q) is neglected. 

Since h g (-,-) is a bilinear form, we can write 

hg((fi + + VO = hg(<f,<p) + (hg(<P,lj)) + hg(lj),tp) + kg {iff , 1^)) 

= hg{<P,V>) + h g - inl (<p,il>). 
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and similarly for /*„(•, •). 

We may now define the multilevel scheme based on the approximation (2.18)- 
(2.23). 

Given U n NM = p n +q n + r", the approximation of ujv.A/(nAf), and an integer 
n u , we define U n w +2 f Uu = p n + 2 '*>‘ +g'*+ 2n “ +r' ,+2n “ by using the following multilevel 
scheme: 


For j = 0, 1, . . . , n u — 1, 
k 2 (l + uk 2 At)(p” +2j+1 ,4u) 


d 2 


-(l + uk 2 At)(-^ P : +2j+ ',0/) u 


a 4 


+ .,At(£ 7 iC +2i+1 ,ifri) =((^-7hK 


% 


a 2 


.n + 2 j 




, </’/ 


+ A<(/j t ,(p n+2 - , ,p n+2j ,V’/)u; + A t(h v ,i»t(p n ,q n + r n ),il'i)«, 

for / = 0, . . . , M p — 4, 


n + 2 j + 1 _ n+2> 

Vv _ ’ 

n+2j+ 1 — n +2j — n. 

1 v ~ v “ 1 r » 

(1 + ^ 2 A7)(^ +22+ \«hL - *A*(|^ 9 " +22+ \ <?,),, = (P g n+2J >')- 

+ Af(/i J (p n+2j ,p , * +2i ), ^i)w + A t(h g , int (p",(q + r) n ), 

for / = 0, . . . , M p — 2, 


n+2j+l _ „«+2.7 

n+2j+l _ r «+2j __ n. 
' 0 9 9 ’ 


(2.24) 


(2.25) 


Q 2 

7k 2 (l + i/fc 2 A<)((Pi> + <Z„)" +2j+2 ,V’<)u> - (1 + vk 2 At)(~{p v + q v ) n+2}+2 , 4uh 

+ uAt + <iv) n+2j+2 ,*!>i S J = - fy2^ Pv + Qv) n+2j+ \ t i’i S j 

+ At(M(P + g) B+2j+1 ,(p + g)" +22+1 ),V-,L 
+ A t{h v , int (p n + q n ,r n ),xl>i) u> , 

for 7 = 0,... , M q — 4, 

n + 2j+2 _ „n + 2ji+l ». 

1 V V V ' 

(2.26) 

(1 + vk 2 At){[pg + 9 g ) n+2j+2 ,<7>/)u; - vAt(-^(pg + 9 9 ) B+2j+2 .^)u- 

= ((p 9 + g 9 r +2i+1 ^/)- 
+ At(ft,((p + «)" +2>+1 ,p n+2i+1 ),^L (2.27) 
+ A t(h.gj nt (p n + q n ,r n ),4>i) u , 

for 7 = 0,... , M q — 2. 

n + 2j + 2 _ r n+2] + \ n 

1 9 3 9 ' 
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Once we obtain p n + 2n “ and g ri+2,,u from above, we compute r n+2 ' lu as follows 


* 2 (1 + 2n u uk' 2 At )(/ , " +2 " u , V-’/ )ai - (1 + 2n„i4--A0(^r" +2, '% V>/)u 




dj/ 2 


+2»,A*(M(P + g)" +2,, “, (p + 9 ) n+2 "“ ),*'•/)„ 


a 2 


n + 2n u 


/")-n) 

2 


(*' 


for / = AT — 3. . . . , M — 4, 


(2.28) 


and 


(1 + 2n „ vk 2 At )( r " +2n “ , *,) w - 2n,,i/A<( T/^ 2 "" , 0,L 

= <r? +2 "',^) w 

( 900 ) 

+ 2n u At(h g ((p + 9)' ,+2 " u , (p + g) n+2 ''u ), ' 

- ((< +2 "“ - «?),*/)„, - 2n.fr 2 Af(9 9 " +2 "«,0,L 

for l = M q - 1, ... , M - 2. 

Note that the computation of the right-hand side of (2.24 )-( 2-25 ) (resp. (2.26)- 
(2.27)) requires only fast Chebyshev transforms (FCT) with 0{M p lo(j 2 {M p )) (resp. 
0( M p log 2 ( M p ) ) ) operations in the normal direction. The nonlinear interaction 
terms h v>rnt and h gant are computed once at the time iteration j — n. Hence, dur- 
ing the 2 n u time iterations described above, FCT with 0{Mlog 2 (M)) are required 
only at j = n and j = n + 2n u . Compared to a classical Galerkin (or tau) approxi- 
mation the multilevel scheme proposed here allows to significantly reduce the CPU 
time needed for channel flow simulations. In the case of forced homogeneous turbu- 
lence, savings of the order of 50-70% have been obtained (see Dubois, Jauberteau 
&; Temam 1995b, 1996). 

3. Incremental unknowns in the finite difference case 

The main idea of the multilevel scheme is to treat the large and small scales dif- 
ferently in numerical simulation. Therefore, it is important to have an appropriate 
decomposition of the flow into different length scales. In Section 3.1, we describe a 
procedure to decompose the solution into large and small scales in finite difference 
method. To illustrate the method, we start by applying the IU’s to the Burger’s 
equation. In Section 3.2, we test the method of separating scales using turbulent 
channel flow database. In Section 3.3, we suggest an algorithm to implement, the 1 
scheme for the turbulent channel flow. 
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3.1 Incremental unknowns on Burger’s equation 
In this section, we start with the two-level IU’s, namely, we decompose the solu- 
tion u into 


u = y + z. 


The second-order IU’s in one-dimensional case are defined as in Chen & Temam 
(1991) by 

V2j = «2 j, 

?2j + l = u 2j + l ~ ^(w2> + 2 + U 2 j). 

Multilevel IU’s can be defined recursively. Three- level IU’s will be defined in Section 
3.2. 

Let us consider the Burger’s equation, 


% + " (<u) = = °- 

When the second order central difference scheme is used for the space derivatives 
and the explicit Euler is used for the time advancing, the finite difference scheme 
reads 


,, n + l ,.n i ,, 

* r/ n ^ 2 1 + + 




At ' AAx iy ~'^ 1 ' ' Ax 2 

Writing y and z components separately, one finds that y satisfies 


u 2j+l — z 2j+\ + < ) (y2) + 2 + y2j)i 

i S 

At ' AAx 
v 


+ tL-k-w ,) 2 - K-.n 


- A?^ 2,+ l “2/2j "b u 2j-l] "b X(x2j,t ), 


and z satisfies 


r n + l 
“2j + l 


~ z Ji±i , _1_ 

At + 2At lvy2jl+2 


+ 


1 

4Ax 


[(y 2 n+1 

[(y”i+2) 2 - (V2j) 2 \ = 


+ y2; +1 )-(y2>+2 + y"j)] 

£^2 [— ^ Z 2j+l] + X( x 2j+l , t ). 


Instead of evaluating z at each time step, we propose to fix z for m steps and then 
evaluate once to save CPU time and memory. Therefore, as m increases, so does the 
saving of CPU time. On the other hand, we are also at the risk of losing accuracy 
as m increases. It is clear that when m = 0, the scheme is the same as the original 
standard method with the fine mesh, while if we never update z and let it be 0, the 
scheme is simply the original standard method in the coarse mesh and u — y. 
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To illustrate how much savings one could obtain by freezing z systematically, we 
list in Table 1 the ratio of the work with freezing c rn times vs. 0 times, with 
the assumption that the work per time step per grid point is independent of the 
mesh size. As an example, if one freezes z for one time step in a three-dimensional 
problem, the work by using IU’s is only 56.25% of that by using the standard finite 
difference method. 

Table 1. Ratio of the work with freezing 2 m times vs. 0 times 


m 

0 

1 

2 

3 


N, 

l-D 

1 

0.75 

0.67 

0.625 


0.50 

2-D 

1 

0.625 

0.5 

0.44 


0.25 

3-D 

1 

0.5625 

0.42 

0.34 


0.125 


We now test this scheme on a model problem, in which we try to recover the 
steady solution u s (x) = /(20a*) — /( 0) + (/( 0) — /(20))a of Burger’s equation, where 

15 

/(*) = exp(cos( k\/k(2.5 + 0.5£)7r/10) — 0.3 sin(0.8ky/ktn/l0)). 

*=i 

The forcing function A *(x,t) — X{x) is calculated by substituting a 9 (x) into the 
equation. Initial condition is taken as u(x,0) = sin(2x) with the boundary con- 
ditions u(0,/) — u(l,t) — 0. By comparing the graphs of u s (a) with N x = 512, 
N x — 256, and N x — 128, one finds that N x = 256 is approximately the minimum 
number of grid points required to adequately resolve u a (x). 

The numerical results using the original scheme and the proposed scheme with 
different m are compared (Fig. 1). For m — 1 to 4, the results are almost iden- 
tical. However, for m — 5, the approximate solution is significantly less accurate. 
Therefore, the proposed scheme has to be used with caution and m can not be too 
large. 

3.2 Small scales in IV 

In the multilevel scheme given in Section 2, a spectral method is used to de- 
compose scales. However, it is not easy to define small scales in finite difference 
methods. In Section 3.1, the small scale component of the flow is defined in the 
context of IU’s. In this section, we examine this concept. 

For simplicity, we will only treat the three-level IU. As is done in Section 2, the 
method starts with separating the length scales of the flow into 

w = / + ff + *\ (3.1) 

where u is the velocity in the streamwise direction, /, g, and r are respectively the 
large, intermediate, and small scales. The definitions of /, g , and r are given below 
(see Fig. 2): 
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FIGURE 1. An exact steady state solution of the Burger’s equation. N=256 is the 
minimum for resolution. 



FIGURE 2. Schematic diagram of separating scales in a finite difference method. 


fa — ^4 < 

1 

</4i+2 = ^4t+2 — ~(u 4l ;+4 + W 4j ) (3.2) 

r 4i + l = **4a + l ~ “(U4i+ 2 + W 4 i) 7 

where i is the index for the streamwise (or wall-normal, or spanwise) direction 
( i = 0, 1 , 2 , N x j4). The wall-normal and spanwise velocities can be defined in a 
similar way. We require the condition 

l/l > \g\ > M (3.3) 

in order to validate the assumption of separating length scales. 

In the present study, the magnitudes of /, g, and r are estimated using the 
database of turbulent channel flow. Turbulent flow in a channel is simulated using 
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DNS. The convection and diffusion terms are integrated in time using a third-order 
Runge-Kutta method and the Crank-Nicolson method, respectively. A second-order 
central difference is used in space. A fractional step method is used to decouple the 
pressure from the velocity. The Reynolds number used is Re T = u T 6/u = 180 and 
the computational domain is 4? tS ( x ) x 26 ( y ) x 4 tt/3 6 (z), where u T is the wall 
shear velocity, 6 is the channel half width, and v is the kinematic viscosity. The 
number of grid points used is 128 (a-) x 129 (j/) x 128 (z). 

Figure 3 shows the energy spectra of the velocity components in the streamwise 
and spanwise directions, where E(fi), E(gj) and E(r t ) are shown at three (/-locations 
(y = 6 , 33 , 177). It is clear that r,’s have the energy of small scales, while g,'s have 
the energy of intermediate scales. Both y, and r, have orders of magnitude smaller 
energies in small wavenumbers as compared to /,. Therefore, the IU’s defined in 
(3.2) properly describe the small and intermediate scales of the velocity. 

3.3 Implementation of IU in turbulent channel flow 

Implementation of IU for the Navier-Stokes equations is very similar to that of 
IU for the Burger s equation (see Section 3.1), once the approximating factorization 
scheme is used (see below). The only difference is the coupling between the velocity 
and the pressure. 

The governing equations for an incompressible flow are 


duj f d dp 1 d du, 

dt dxj U,U} ~ dx, + Redx]dx]' 


(3.4) 


The integration method used to solve (3.4) and (3.5) is based on a semi-implicit 
fractional step method, i.e., third-order Runge-Kutta method for the convection 
terms and Crank-Nicolson method for the diffusion terms: 


u k — u k 1 


—i a k + 0k)Li(u k 1 ) + flkLi(u k — u k *) 
-7i kN i (u k ~ 1 )~CkN,(u k - 2 ), 

2,* _ 1 du k 


v z <j> 


u k — u k 


At dxi ’ 

- d<t>k 
dx , 7 


(3.6) 


(3.7) 


At dx , 7 ( 3 - 8 ) 

where Li and X, are the diffusion and convection terms of (3.4), k = 1,2,3, and 


a,=A = -, 

<*2 = 02 = — , 
15 

, 1 

«3 = &3 = ", 

0 


8 

71 = 15 7 Cl= ° 


72 12 7 C ' 2 60 


17 
30 

3 . 5 

73 = 4 7 = 
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Here, ( a k + (h )p k = 4> k - (Af/?*/ Re)V 2 <j> k . 

Rearranging (3.6) in delta form (6u k = uf — u k ~ l ) gives 

(1 - Atfa±V 2 )6u k = At [(a* + flkJXifti*- 1 ) - itJViiv*- 1 ) - ( k N,(u k ~ 2 )}. 
Approximating factorization of this equation gives 

(1 ~ ~ ^‘‘Redf^ 1 ~ ^ k Red^^ 

= Af[(a t + ih)Li(u k ~ l ) - lk Ni{u k - 1 ) - C**Vj(n‘“ 2 )] (3,9) 

= Ri{u k -' ,u k ~ 2 ). 

Let us define \, as 

( 3 . 10 ) 

Then, (3.9) becomes 

(1 - Af/3 fc — — = -R;. (3.11) 

For simplicity, we only focus on the velocity in the streamwise component. Note 
that in turbulent channel flow the periodic boundary conditions are applied in the 
streamwise and spanwise directions (,r, z) and the no-slip condition is applied in the 
wall- normal direction (y). 

Now, let us decompose \ (streamwise component of \,) into three different scales 
as was introduced in Section 3.2: 


x= f + g + r 


(3.12) 


As a first step, (3.11) is approximated at each fourth grid point using a second- 
order central difference scheme: 


\4< — r(\4i+l — 2 x4i + \4*-l ) = Rl 4i , (3.13) 

where T = Re Aa: 2 ). 

Using a similar relation to (3.2), it can be easily shown that (3.13) becomes 


— “\4t+4 + (1 + ^)\4i — j\4i-4 

_ _,1 1 

— Hi 4, + F(-y 4l + 2 + ^#4i-2 + r 4r+l + ^4i-l )• 


(3.14) 


The \ at every fourth grid point is obtained by solving (3.14). The \4t±i and \' 4 t ±2 
are updated with the newly obtained \4i from (3.14): e.g., 

\4i+2 = </li+2 + y(X4i T \4i+4 ) 

\4i+l — r 4*-hl + ^(\4i + \4f+2), 


(3.15) 
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where g and r are frozen for the periods of A t g and A £ r , respectively. A t g and At r 
are called as the frozen times for the intermediate and small scales, respectively. 

As a second step, (3.11) is approximated at every other grid point at t = lAt g (/ 
is an integer) using a second-order central difference: 

r r 

- 2 X 2 i +2 + (1 + r)\ 2 i - 2 * 2t - 2 = R-U i + r(r 2 ,-+i + r 2 i- 1). (3.16) 


The x at every other grid point is obtained by solving (3.16). The X 2 i±i are then 
updated with the frozen r, and gn± 2 are updated: e.g., 


1, 

X2i+1 = r 2i+l + 2^ 2 ' + X2i+2) 

(3.17) 

9ii+2 = X4i+2 ~ "*■ X4i+4)- 

(3.18) 

As a third step, (3.11) is approximated at every point at t = lAt r : 


—■ TXi+i + (1 + 2r)xi — r^i-i = Ri r 

(3.19) 

The \ at every grid point is obtained by solving (3.19). The r 2 ,±i are 

then updated 

as 


r 2»+l = X2i+1 — «(X2i + X2»+2 )* 

(3.20) 


Once x’s are obtained at either 4 i, 2i, or i points, similar procedures are applied 
to the other two directions. It is straightforward to extend the procedure described 
above in the spanwise and wall-normal directions. At the end of these procedures, 
the streamwise velocity is obtained. Again, the same procedure can be easily applied 
to the other two velocity components. 

Let us write the numerical algorithm of IU: 

1. Start with an initial velocity field u° or a previous time step u n, * _1 = u”' 1 . 

2. Solve the discretized momentum equations at (4i,4j,4fc) grid points (similar to 
(3.14)) to obtain u at (4i,4j,4fc) points. 

3. Update u at non-(4f, 4j, 4k) points with frozen g and r (see (3.15)). 

4. If t = /Af 5 , go to Step 5. If not, go to Step 2. 

5. Solve the discretized momentum equations at (2i,2j,2Ar) grid points (similar to 
(3.16)) to obtain u at (2i,2j, 2k) points. 

6. Update u at (2 i ± 1,2 j ± 1,2 k ± 1) points with frozen r and also update g at 
(4t ±2,4j ±2,4* ±2) points (see (3.17) - (3.18)). 

7. If t = /At r , go to Step 8. If not, go to Step 2. 

8. Solve the discretized momentum equations at all the grid points (similar to (3.19)) 
to obtain u at all points. 

9. Update r at (2 i ± 1,2 j ± l,2Ar ± 1) points (see (3.20)). 

10. Solve the Poisson Eq. (3.7) at all points, update the velocity (3.8), and go to Step 

2 . 
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Note that it is not necessary for us to decompose the velocity into the same 
levels of scales in all the directions. That is, one may decompose the flow into two 
scales in the wall-normal direction and three scales in the streamwise and spanwise 
directions. 

The interpolation used in obtaining the neighboring velocity (e.g., (3.15)) dete- 
riorates the momentum conservation property, and the mass conservation is easily 
violated unless the Poisson Eq. (3.7) is solved at each time step. However, the re- 
quirement of the Poisson solution at each time step clearly diminishes the advantage 
of using the IU method. 

The modification and application of the present multilevel scheme to the turbulent 
channel flow are in progress and will be reported in the future. 
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